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Abstract 

Recovering the shape of an object from two views fails at occluding contours of smooth objects 
because the extremal contours are view dependent. For three or more views, shape recovery is 
possible, and several algorithms have recently been developed for this purpose. We present a new 
approach to the multiframe stereo problem which does not depend on differential measurements in 
the image, which may be noise sensitive. Instead, we use a linear smoother to optimally combine 
all of the measurements available at the contours (and other edges) in all of the images. This allows 
us to extract a robust and dense estimate of surface shape, and to integrate shape information from 
both surface markings and occluding contours. 
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1 Introduction 

Most visually-guided systems require representations of surfaces in the environment in order to 
integrate sensing, planning, and action. The task considered in this paper is the recovery of the 3D 
structure (shape) of objects with piecewise- smooth surfaces from a sequence of profiles taken with 
known camera motion. The profile (also known as the extremal boundary or occluding contour) is 
defined as the image of the critical set of the projection map from the surface to the image plane. 
Since profiles are general curves in the plane without distinguished points, there is no a priori 
pointwise correspondence between these curves in different views. However, given the camera 
motion, there is a correspondence based on the epipolar constraint. For two images, i.e., classical 
binocular stereo, this epipolar constraint is a set of straight lines which are the intersection of the 
epipolar planes with the image plane. The epipolar plane through a point is determined by the 
view direction at that point and the instantaneous camera translation direction. 

In the case of contours that are not view dependent, e.g., creases (tangent discontinuities) and 
surface markings, many techniques have been developed for recovering the 3D contour locations 
from two or more images under known camera motion [Marr and Poggio, 1979; May hew and 
Frisby, 1980; Arnold, 1983; Bolles et al, 1987; Baker and Bolles, 1989; Matthies et al, 1989]. 
Techniques have also been developed for simultaneously estimating contour locations and camera 
positions [Tsai and Huang, 1984; Faugeras et al, 1987; Horn, 1990]. However, for smooth curved 
surfaces, the critical set which generates the profile is different for each view. Thus, the triangulation 
applied in two-frame stereo will not be correct along the occluding contour for smooth surfaces. 
For the same reason, it is often not possible to determine the camera motion from the images unless 
some assumptions are made either about the surface or the motion [Arborgast and Mohr, 1992; 
Giblin et al, 1992]. On the other hand, the fact that the critical sets sweep out an area means that 
the connectivity of the surface points can be determined, i.e., one obtains a surface patch rather 
than a set of points. 

The problem of reconstructing a smooth surface from its profiles has been explored for known 
planar motion by Giblin and Weiss [Giblin and Weiss, 1987] and subsequently for more general 
known motion by Vaillant and Faugeras [Vaillant, 1990; Vaillant and Faugeras, 1992] and Cipolla 
and Blake [Blake and Cipolla, 1990; Cipolla and Blake, 1990; Cipolla and Blake, 1992]. These 
approaches are either based on a differential formulation and analysis, or they use curve fitting but 
still only use three frames. First order temporal derivatives are usually computed as differences from 
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pairs of frames, and second order derivatives from triples. Unfortunately, determining differential 
quantities reliably from real images in this way is difficult. Even fitting curves to data from three 
images can be unsatisfactory. This has led Cipolla and Blake to use relative measurements in 
order to cancel some of the error due to inadvertent camera rotation. Their approach approximated 
image contours with B-snakes which require initialization for each contour that is tracked. In 
addition, B-snakes implicitly smooth the contours in the image. Since the recovery of 3D points is 
a linear problem (as we will show in this paper), the smoothing can be done in 3D on the surface, 
where more context can be used in the detection of discontinuities so that detailed structure can be 
preserved. 

It is natural to consider surface reconstruction as an optimal estimation problem. To overcome 
the limitations of previous algorithms, the approach we develop in this paper applies standard 
techniques from estimation theory (Kalman filtering and smoothing) to make optimal use of each 
measurement without computing differential quantities. First, we derive a linear set of equations 
between the unknown shape (surface point positions and radii of curvature) and the measurements. 
We then develop a robust linear smoother ([Gelb, 1974; Bierman, 1977]) to compute statistically 
optimal current and past estimates from the set of contours. Smoothing allows us to combine 
measurements on both sides of each surface point. 

Our technique produces a complete surface description, i.e., a network of linked 3D surface 
points, which provides us with a much richer description than just a set of 3D curves. Some parts 
of the surface may never appear on the profile. In some cases this is due to occlusion either by the 
same surface or another one. In other cases, it is due to limitations of the camera trajectory [Giblin 
and Weiss, 1994]. Since the method presented here also works for arbitrary surface markings 
and creases, a larger part of the surface can be reconstructed than from occluding contours of the 
smooth pieces alone. Our approach also addresses the difficult problem of contours that merge and 
split in the image, which must be resolved if an accurate and complete 3D surface model is to be 
constructed. 

The method we develop has applications in many areas of computer vision, computer aided 
design, and visual communications. The most traditional application of visually based shape 
recovery is in the reconstruction of a mobile robot's environment, which allows it to perform 
obstacle avoidance and planning tasks [Curwen et ah, 1992]. Visually based shape recovery can 
also be used to develop strategies for robotics grasping and manipulation tasks, or as an off- 
line technique to automatically "learn" object descriptions for object or pose recognition tasks. 
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In less traditional applications, our system could be used to perform reverse engineering, i.e., 
to automatically acquire 3D computer aided design (CAD) description of real-world objects or 
prototypes, or even to construct a "3D fax" for transmitting 3D object descriptions and images 
between desktops. 

Our paper is structured as follows. We begin in Section 2 with a description of our edge 
detection, contour linking, and edge tracking algorithms. In Section 3, we discuss the estimation of 
the epipolar plane for a sequence of three or more views. Section 4 presents the linear measurement 
equations which relate the edge positions in each image to the parameters of the circular arc 
being fitted at each surface point. Section 5 then reviews robust least squares techniques for 
recovering the shape parameters and discusses their statistical interpretation. Section 6 shows 
how to extend least squares to a time-evolving system using the Kalman filter, and develops the 
requisite forward mapping (surface point evolution) equations. Section 7 extends the Kalman filter 
to the linear smoother, which optimally refines and updates previous surface point estimates from 
new measurements. Section 9 presents a series of experiments performed both on noisy synthetic 
contour sequences and on real video images. We close with a discussion of the performance of our 
new technique and a discussion of future work. 

2 Contour detection and tracking 

The problem of edge detection has been extensively studied in computer vision [Marr and Hildreth, 
1980; Canny, 1986]. The choice of edge detector is not crucial in our application, since we 
are interested mostly in detecting strong edges such as occluding contours and visible surface 
markings. 1 For our system, we have chosen the steerable filters developed by Freeman and 
Adelson [Freeman and Adelson, 1991], since they provide good angular resolution at moderate 
computation cost, and since they can find both step and peak edges. We have used both the Gi 
and (G2, H2) sets of filters, with the default parameters suggested by Freeman and Adelson. An 
example of our edge detector operating on the input image in Figure la is shown in Figure lb. 

Once discrete edgels have been detected, we use local search to link the edgels into contours. 
We find the two neighbors of each edgel based on proximity and continuity of orientation. Note 
that in contrast to some of the previous work in reconstruction from occluding contours [Cipolla 

'Unlike many edge detection applications, however, our system provides us with a quantitative way to measure the 
performance of an edge detector, since we can in many cases measure the accuracy of our final 3D reconstruction. 
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(c) (d) 

Figure 1: Input processing steps 
(a) sample input image (dodecahedral puzzle), (b) estimated edgels and orientations (maxima in 
| Gi | 2 ), (c) tracked edgels, (d) correspondence of points on the occluding contours using the epipolar 
constraint. 
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and Blake, 1990; Cipolla and Blake, 1992; Blake et ai, 1993], we do not fit a smooth parametric 
curve to the contour since we wish to directly use all of the edgels in the shape reconstruction, 
without losing detail. 2 The curve fitting problem is essentially one of detecting outliers. Since the 
3D reconstruction provides more context, smoothing in 3D should be preferred. 

We then use the known epipolar line constraints (Section 3) to find the best matching edgel in 
the next frame. Our technique compares all candidate edgels within the epipolar line search range 
(defined by the expected minimum and maximum depths), and selects the one which matches most 
closely in orientation, contrast, and intensity (see Figure lc). Once an initial estimate for the 3D 
location of an edgel has been computed, the search range can be dramatically reduced (see Section 
5.3). 

Since contours are maintained as a list of discrete points, it is necessary to resample the edge 
points in order to enforce the epipolar constraint on each track. We occasionally start new tracks if 
there is a sufficiently large (2 pixel wide) gap between successive samples on the contour. While we 
do not operate directly on the spatiotemporal volume, our tracking and contour linking processes 
form a virtual surface similar to the weaving wall technique of Harlyn Baker [Baker, 1989] . Unlike 
Baker's technique, however, we do not assume a regular and dense sampling in time. 

3 Reconstructing surface patches 

The surface being reconstructed from a moving camera can be parametrized in a natural way by two 
families of curves [Giblin and Weiss, 1987; Cipolla and Blake, 1990]: one family consists of the 
critical sets on the surface; the other is tangent to the family of rays from the camera focal points. 
The latter curves are called epipolar curves and together with the critical sets form the epipolar 
parametrization. This parametrization can always be used except when the profile is singular or 
when the normal to the surface is perpendicular to the direction of camera translation [Giblin and 
Weiss, 1994]. For a pair of stereo images, each viewing direction together with the translation 
vector from one camera center to the other determines a plane called an epipolar plane. The same 
construction holds in the case of motion: in the limit, as the time between samples goes to zero, 
the plane determined by a view direction and the camera translation velocity will also be called an 
epipolar plane. For a more detailed discussion of epipolar curves see [Giblin and Weiss, 1994]. 

2 However, we do perform a small amount of curvature-dependent smoothing along the curves to reduce noise. This 
can be viewed as part of the edge extraction stage. 
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The problem is that any smooth surface reconstruction algorithm which is more than a first order 
approximation requires at least three images and, that in general, the three corresponding tangent 
rays will not be coplanar. However, there are many cases when this will be a good approximation. 
One such case is when the camera trajectory is almost linear. If the camera trajectory is linear, then 
the epipolar planes form a pencil of planes containing that line. Under orthographic projection, if 
the camera motion is planar, then all of the epipolar curves will be planar as well. 

Cipolla and Blake [Cipolla and Blake, 1990; Cipolla and Blake, 1992] and Vaillant and Faugeras 
[Vaillant, 1990; Vaillant and Faugeras, 1992] noticed that to compute the curvature of a planar curve 
from three tangent rays, one can determine a circle which is tangent to these rays. See Figure 2. 
The assumption that one needs to make is that the surface remains on the same side of the tangent 
rays. This is true for intervals of the curve which do not have a singularity or zero curvature. 

Given three or more edgels tracked with our technique, we would like to compute the location 
of the surface and its curvature by fitting a circular arc to the lines defined by the view directions 
at those edgels. In general, a non-singular space curve will have a unique circle which is closest 
to the curve at any given point. This is called the osculating circle, and the plane of this circle 
is called the osculating plane. It is easy to see that if the epipolar curve is non-singular, then the 
epipolar plane is an estimate of its osculating plane [Cipolla and Blake, 1992], and the lines defined 
by the view directions are close to this plane and can be projected onto it. The accuracy of the 
computation of the radius of curvature depends on the conditioning of this projection. Since in the 
limit, the epipolar plane is the osculating plane for the epipolar curve, the epipolar curves should 
be the most robust to reconstruct by projecting onto this plane. 

The relationship between the curvature of a curve such as the epipolar curve and the curvature 
of the surface is determined by the angle between the normal to the curve and the normal to the 
surface. The curvature of the curve scaled by the cosine of this angle is the normal curvature. The 
curvature of a surface can be thought of as a function which assigns to every tangent direction v 
a value which is the curvature k v of the normal slice in that direction. If v is the tangent to the 
epipolar curve, k epi is the curvature of the epipolar curve, and 4> is the angle between the epipolar 
plane and the plane containing v and the surface normal, then the relationship among them is given 
by the equation 

k v = k epi cos^ (1) 

This gives Meusnier's Theorem which says that the normal curvature is the same for all curves on 
the surface with a given tangent direction. Since the normal to the surface can be determined from 
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c 




Figure 2: Local coordinate axes and circle center point calculation 
In this figure, the points p ; and q ; are coincident. In general, the q ; will lie somewhere along the 
tangent rays, and the p ; will be the points of tangency to the osculating circle. 



the image, the normal curvature can be obtained from the epipolar curve. 

4 Measurement equations 

Once we have selected the epipolar plane as the reconstruction plane for fitting the circular arc, we 
must compute the set of lines in this plane which should be tangent to the circle. This can be done 
either by projecting the 3D lines corresponding to the linked edgels directly onto the plane, or by 
intersecting the tangent planes (defined by the edgels and their orientations) with the reconstruction 
plane. 

We represent the 3D line corresponding to an edgel in frame i by a 3D point q; and a direction t; . 
The point q; is chosen to be the intersection of the viewing ray with a reference plane z = z$. The 
direction is given by t; = A/"(q; — c;), where c; is the camera center and Af() normalizes a vector. 
We choose one of these lines as the reference frame (no, to) centered at q 0 (where iio = t 0 x n ep i), 
e.g., by selecting the middle of n frames for a batch fit, or the last frame for a Kalman filter. This 
line lies in the reconstruction plane defined by n epi . 
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If we parameterize the osculating circle by its center c = (x c ,y c ) and radius r (Figure 2), we 
find that the tangency condition between line i and the circle can be written as 

c { x c + s { y c + r = di (2) 

where c; = t; • to, = — t; • iio, and di = (q; — q 0 ) • n;. Thus, we have a linear estimation problem 
in the quantities (x c ,y c ,r) given the known measurements (q, s;, c?;). This linearity is central to 
the further developments in the paper, including the least squares fitting, Kalman filter, and linear 
smoother, which we develop in the next three sections. 



5 Least squares fitting 



While in theory the equation of the osculating circle can be recovered given the projection of three 
non-parallel tangent lines onto the epipolar plane, a much more reliable estimate can be obtained 
by using more views. Given the set of equations (2), how can we recover the best estimate for 
(x c ,y c ,r)7 Regression theory [Albert, 1972; Bierman, 1977] tells us that the minimum least 
squared error estimate of the system of equations Ax = d can be found by minimizing 



e = | Ax - d| 2 = Y,( a i ' x _ d i) 2 - 

i 

This minimum can be found by solving the set of normal equations^ 

(A T A)x = A T d 



(3) 



(4) 



or 



A statistical justification for using least squares will be presented shortly (Section 5.1). 
In our circle fitting case, a ; = (c ; , s;, 1), x = (x c , y c , r), and the normal equations are 



Ei Ci Ei Ei Ci 

Ei Si^i Ei &i Ei &i 
Ei c i Ei s i Ei 1 



Xq 




Ei Cidi 


Vc 




Ei &idi 


r 




Ei^i 



(5) 



Alternative techniques for solving the least squares problem include singular value decomposition [Press et ah, 
1986] and Householder transforms [Bierman, 1977]. 
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If we solve the above set of equations directly, the estimates for x c and r will be very highly 
correlated and both will be highly unreliable (assuming the range of viewpoints is not very large). 
This can be seen both by examining Figure 2, where we see that the location of c is highly sensitive 
to the exact values of the ck, or by computing the covariance matrix P = (A T A) _1 (Section 5.1). 

We cannot do much to improve the estimate of r short of using more frames or a larger camera 
displacement, but we can greatly increase the reliability of our shape estimate by directly solving 
for the surface point (x s , y s ), where x s = x c + r and y s = y c . 4 The new set of equations is thus 

CiX s + s { y s + (1 - Ci)r = d { . (6) 

While there is still some correlation between x s and r, the estimate for x s is much more reliable 
(Section 5.1). Once we have estimated (x s ,y s ,r), we can convert this estimate back to a 3D surface 
point, 

Po = qo + x s h 0 + y s t 0 , (7) 

a 3D center point 

c = q 0 + (a;, - r)n 0 + y s t 0 = Po - rh 0 , (8) 
or a surface point in the i th frame 

Pi = c + rhi = po + r(hi - n 0 ), (9) 

where 

fi» = x n e pi 

is the osculating circle normal direction perpendicular to line U (Figure 2). 
5.1 Statistical interpretation 

The least squares estimate is also the minimum variance and maximum likelihood estimate (optimal 
statistical estimate) under the assumption that each measurement is contaminated with additive 
Gaussian noise [Bierman, 1977]. If each measurement has a different variance oj, we must weight 
each term in the squared error measure (3) by W{ = a[ 2 , or, equivalently, multiply each equation 
a; • x = di by a^ 1 . 

4 While the point (x„, y s ) will not in general lie on the line (q 0 , to), the tangent to the circle at (x„, y s ) will be 
parallel to to- 
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In our application, the variance of c? ; , of, can be determined by analyzing the edge detector 
output and computing the angle between the edge orientation and the epipolar line 

o 2 = a 2 J(h ■ n epi ) 2 = o e 2 /(l - (m, • n epi ) 2 ), 

where o e is the variance of q; along the surface normal ihj. This statistical model makes sense if 
the measurements d{ are noisy and the other parameters (c;, s;) are noise-free. This is a reasonable 
assumption in our case, since the camera positions are known but the edgel locations are noisy. 
The generalization to uncertain camera locations is left to future work. 

When using least squares, the covariance matrix of the estimate can be computed from P = 
(A T A) _1 . We can perform a simple analysis of the expected covariances for n measurements 
spaced 6 apart. Using Taylor series expansions for q = cos i6 and s; = s'mi6, and assuming that 



. m], n 


= 2m + 1, we obtain the covariance matrices 








60" 4 0 -66~ 4 ' 
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where P 3 C is the 3 point covariance for the center-point formulation, and P* is the 3 point covariance 
for the surface-point formulation. As we can see, variance of the surface point local x estimate 
is four orders of magnitude smaller than that of the center point. Similar results hold for the 
overdetermined case (n > 3). Extending the analysis to the asymmetrical case, i £ [0 . . . 2m], we 
observe that the variance of the x s and y s estimates increases. 

5.2 Robustifying the estimate 

To further improve the quality and reliability of our estimates, we can apply robust statistics 
to reduce the effects of outliers which are due to grossly erroneous measurements as well as 
large changes in the surface curvature [Huber, 1981]. Many robust techniques are based on first 
computing residuals, r ; = d{ — a ; • x, and then re-weighting the data by a monotonic function 

W)- 2 = *rVkil) 

or throwing out measurements whose |r;| ^> o;. Alternatively, least median squares can also be 
used to compute a robust estimate, but at an increased complexity. 
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In our application, outliers occur mainly from gross errors in edge detection (e.g., when adjacent 
edges interfere) and from errors in tracking. Currently, we compute residuals after each batch fit, 
and keep only those measurements whose residuals fall below a fixed threshold. 

5.3 Predicting 2D locations for tracking 

Once a 3-D estimate for an edgel location has been computed, this can be used to predict where the 
edgel would appear in the next frame, and hence to improve the correspondence produced by the 
tracking stage. When no 3-D information is available, we project the viewing ray passing through 
a 2-D edgel into the next frame to give us the epipolar search line. We use the intersection of the 
viewing ray with a minimum and maximum depth plane to determine the endpoints which limit the 
search range. 

When a 3-D position estimate is available, we project the 3-D position and covariance estimate 
into the new reconstruction plane. The position on the screen of the edgel then gives us the middle 
of the search range, while a multiple of the standard deviation in the local x direction (which is 
parallel to the image plane and in the reconstruction plane and hence along the epipolar line) times 
the epipolar line determines the limits of the search range. More formally, the endpoints of the 
search line are 

Vi(pi± a<r x h 0 ) 

where V{ projects points in 3-D onto the ith frame, and al is the variance in the local x direction. 

Our approach is similar in spirit to the validation gate approach used by Blake et al. for 
their Kalman-filter snake tracking [Blake et al, 1993]. Even more sophisticated data association 
techniques could be used to disambiguate multiple intersecting tracks [Bar-Shalom and Fortmann, 
1988]. 

6 Kalman filter 

The Kalman filter is a powerful technique for efficiently computing statistically optimal estimates of 
time- varying processes from series of noisy measurements [Gelb, 1974; Bierman, 1977; Maybeck, 
1979]. In computer vision, the Kalman filter has been applied to diverse problems such as motion 
recovery [Rives et al., 1986], multiframe stereo [Matthies et al., 1989], and pose recovery [Lowe, 
1991]. In this section, we develop a Kalman filter for contour-based shape recovery in two parts: 
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first, we show how to perform the batch fitting of the previous section incrementally; second, we 
show how surface point estimates can be predicted from one frame (and reconstruction plane) to 
another. 

The update or data processing part of a Kalman filter takes a current estimate 5q with its 
associated covariance P ; and produces an updated estimate x; and covariance P ; by processing a 
single measurement 

di = a; • X;. (10) 
The traditional Kalman filter formulation [Gelb, 1974] first computes a Kalman gain matrix 

K i = P i a i (afP i a i + ^ 2 )- 1 , (11) 

where cr? is the variance associated with measurement i. It then increments the state estimate by 
adding a weighted residual 

Xi = X; + Ki(di - a; • Xi). (12) 

and decrements the covariance matrix 

Pi = Pi - Kiaf Pi. (13) 

Applying this Kalman filter to our circular arc fitting task is straightforward, since each of our 
tangent lines is of the required form (10), = a» • x». More numerically stable or computationally 
efficient forms of the Kalman filter have also been developed [Bierman, 1977], but we have not yet 
implemented them to see if they improve our performance. 

The update part of the Kalman filter is derived directly from the measurement equation (2) 
[Gelb, 1974] . It provides an incremental technique for estimating quantities in a static system, e.g., 
for refining a set of (x s ,y s , r) estimates as more edgels are observed. For our application, however, 
we need to produce a series of surface points which can be linked together into a complete surface 
description. If we were using batch fitting, we would perform a new batch fit centered around each 
new 2D edgel. Instead, we use the complete Kalman filter, since it has a much lower computational 
complexity. The Kalman filter provides a way to deal with dynamic systems where the state x; is 
evolving over time. We identify each measurement Xi with the surface point (x s ,y s ,r) whose local 
coordinate frame is given by (n,, tj, n ep i) centered at qi in frame i. 

The second half of the Kalman filter requires a system model or process model which predicts 
the evolution of state variables over time [Gelb, 1974]. For our smooth surface model, we assume 
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that r (the third component of x) can vary slowly over time, but that the other two components 
have no associated process noise, i.e., s = (0, 0, ^). 

The overall sequence of processing steps is therefore the following. Initially, we perform a 
batch fit to n > 3 frames, using the last frame as the reference frame. Next, we convert the local 
estimate into a global 3D position (7) and save it as part of our final surface model. We use this 3D 
estimate to construct a reduced search range for edgels during the tracking phase. Then, we project 
the 3D surface point and its radius onto the next frame, i.e., into the frame defined by the next 
2D edgel found by the tracker. 5 Then, we update the state estimate using the local line equation 
and the Kalman filter updating equations. We repeat the above process (except for the batch fit) 
so long as a reliable track is maintained (i.e., the residuals are within an acceptable range). If the 
track disappears or a robust fit is not possible, we terminate the recursive processing and wait until 
enough new measurements are available to start a new batch fit. 

7 Linear smoothing 

The Kalman filter is most commonly used in control systems applications, where the current 
estimate is used to determine an optimal control strategy to achieve a desired system behavior 
[Gelb, 1974]. In certain applications, however, we may wish to refine old estimates as new 
information arrives, or, equivalently, to use "future" measurements to compute the best current 
estimate. Our shape recovery application falls into this latter category, since the accuracy of the 
estimate depends on the range of viewing angles for the measurements, and this can be increased 
by taking measurements from both sides of the 3D curve corresponding to a given visible occluding 
contour. In addition, it should be noted that if the curvature of the epipolar curve is not constant, 
then for each interval over which it is monotonic, filtering rather than smoothing will introduce a 
bias. 

The generalization of the Kalman filter to update previous estimates is called the linear smoother 
[Gelb, 1974]. The smoothed estimate of x; based on all the measurements between 0 and N is 
denoted by x^. Three kinds of smoothing are possible [Gelb, 1974]. In fixed-interval smoothing, 
the initial and final times 0 and N are fixed, and the estimate x^ is sought, where i varies from 0 
to N. In this case, each point in the model is estimated from all of the data in a track. Infixed-point 

5 For even higher accuracy, we could use the 2D projection of our 3D surface point as the input to our tracker. 
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smoothing, i is fixed and x;^ is sought as N increases. Each point is updated as new data is 
obtained, i.e. there is a separate smoother for each point. In fixed-lag smoothing, xn-l\n is sought 
as N increases and L is held fixed. This has the advantage that each point is estimated within a 
fixed amount of time from when it appears on the profile, and it is only estimated once. Since the 
lag time is non-zero, information on both sides of the critical set are used. 

For surface shape recovery, both fixed-interval and fixed-lag smoothing are of interest. Fixed- 
interval smoothing is appropriate when shape recovery is performed off-line from a set of prede- 
termined measurements. The results obtained with fixed-interval smoothing should be identical 
to those obtained with a series of batch fits, but at a much lower computational cost. The fixed- 
interval smoother requires a small amount of overhead beyond the regular Kalman filter in order to 
determine the optimal combination between the outputs of a forward and backward Kalman filter 
[Gelb, 1974; Bierman, 1977]. 

For our contour-based shape recovery algorithm, we have developed a new fixed-lag smoother, 
which, while sub-optimal, allows us to predict the position of the contour in successive images and 
simplifies the tracking problem. Our fixed-lag smoother begins by computing a centered batch fit 
to n(> 3) frames. The surface point is then predicted from frame i — 1 to frame i as with the 
Kalman filter, and a new measurement from frame i + L, L = [n/2\ is added to the predicted 
estimate. The addition of measurements ahead of the current estimate is straightforward using the 
projection equations for the least-squared (batch) fitting algorithm. 

Our modified fixed-lag smoother and the optimal fixed-lag smoother incorporate the same 
information into the current estimate, but use slightly different relative weightings of the data. 
Intuitively, the optimal smoother weights the data most heavily towards the middle (inversely 
proportional to the distance from the current estimate), while our modified smoother weights the 
data most heavily towards the front (most recent measurement). For systems where the process 
noise a] is much smaller than the measurement noise a}, the results should be similar. We examine 
the relative performance of the batch estimator, Kalman filter, and sub-optimal linear smoother in 
Section 9. 

8 Building a complete surface description 

The batch fitting, Kalman filter, and linear smoothers all produce a series of surface point estimates, 
one for each input image. Because our reconstruction takes place in object space, features such as 
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surface marking and sharp ridges are stationary in 3D (and have r = 0). For these features, we 
would prefer to produce a single time-invariant estimate. While the detection of stationary features 
could be incorporated into the Kalman filter or smoother itself, we currently defer this decision 
to a post-processing stage, since we expect the estimates of position and radius of curvature to be 
more reliable after the whole sequence has been processed. The post-processing stage collapses 
successive estimates which are near enough in 3D (say, less than the spacing between neighboring 
sample points on the 3D contour). It adjusts the neighbor (contour) and temporal (previous/next) 
pointers to maintain a consistent description of the surface. 

To fit a complete surface to the data while interpolating across small gaps, a variety of techniques 
could be used. Traditionally, physically-based deformable models [Terzopoulos and Metaxas, 
1991; Pentland and Sclaroff, 1991] have been used to fit such sparse and incomplete 3-D data. An 
alternative, which does not suffer from the restrictions on topology imposed by previous techniques, 
is to use a self-triangulating system of particles to model and interpolate the surface [Szeliski etal, 
1993]. We plan to investigate the intergration of this system with our multiframe stereo algorithm 
in future work. 

9 Experimental results 

To determine the performance of our shape reconstruction algorithm, we generated a synthetic 
motion sequence of a truncated ellipsoid rotating about its z axis (Figure 3). The camera is oblique 
(rather than perpendicular) to the rotation axis, so that the epipolar curves are not planar, and the 
reconstruction plane is continuously varying over time. We chose to use a truncated ellipsoid since 
it is easy to analytically compute its projections (which are ellipses, even under perspective), and 
since its radius of curvature is continuously varying (unlike, a cylinder). 

When we run the edge images through our least-squares fitter or Kalman filter/smoother, we 
obtain a series of 3D curves. The curves corresponding to the surface markings and ridges (where 
the ellipsoid is truncated) should be stationary and have 0 radius, while the curves corresponding 
to the occluding contour should continuously sweep over the surface. 

We can observe this behavior using a three-dimensional graphics program we have developed 
for displaying the reconstructed geometry. This program allows us to view a series of reconstructed 
curves either sequentially (as an animation) or concurrently (overlayed in different colors), and to 
vary the 3D viewing parameters either interactively or as a function of the original camera position 
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Figure 3: Four images from synthetic truncated ellipsoid sequence. 
The top and left hand side are truncated (cut off), while the front and back sides are inscribed with 
an ellipse (surface marking). 



Figure 4: Top view of reconstructed 3D curves. 
The surface markings and ridges are stationary, while the occluding contours (ellipse) sweeps 
around the object. 



for each frame. Figure 5 shows all of the 3D curves overlayed in a single image. As we can see, 
the 3D surface is reconstructed quite well. The left hand pair of images shows an oblique and top 
view of a noise-free data set, using the linear smoother with n = 7 window size. The right-hand 
pair shows a portion of the reconstructed surface, showing both the profile and epipolar curves. 

To obtain a quantitative measure of the reconstruction algorithm performance, we can compute 
the root median square error between the reconstructed 3D coordinates and the true 3D coordinates 
(which are known to the synthetic sequence generating program). Table 1 shows the reconstruction 
error and percentage of surface points reconstructed as a function of algorithm choice and various 
parameter settings. The table compares the performance of a regular 3-point fit with a 7-point 
moving window (batch) fit, and a linear fixed-lag smoother with n = 7. Results are given for the 
noise-free and cr ; = 0. 1 pixels case. The different columns show how by being more selective about 
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Figure 5: Oblique and top view of reconstructed 3D surface (all 3D curves superimposed). 
The left pair shows only the reconstructed profile curves, while the right pair shows the profiles 
linked by the epipolar curves (only a portion of the complete meshed surface is shown for clarity). 
A total of 72 images spaced 5° apart were used. 

which 3D estimates are considered valid (either by requiring more frames to have been successfully 
fit, or lowering the threshold on maximum covariance), a more reliable estimate can be obtained at 
the expense of fewer recovered points. For noisefree data, the 3 point algorithm is better because it 
is less sensitive to curvature variation. However, for noisy data, the 7 point algorithms are better, 
with batch fitting performing slightly better than linear smoothing. 

We have also applied our algorithm to the four real image sequences shown in Figure 6. 
These sequences were obtained by placing an object on a rotating mechanized turntable whose 
edge has a Gray code strip used for reading back the rotation angle [Szeliski, 1991; Szeliski, 
1993]. The camera motion parameters for these sequences were obtained by first calibrating the 
camera intrinsic parameters and extrinsic parameters to the turntable top center, and then using the 
computed turntable rotation. Figure 7 shows the edges extracted from each of these images. 

Figure 8 shows two views of each set of reconstructed 3D curves. We can see that the overall 
shape of the objects has been reconstructed quite well. We show only the profile curves, since 
the epipolar curves would make the line drawing too dense for viewing at this resolution. Figure 
9 shows both the profile curves and the epipolar curves for selected portions of the soda can and 
coffee objects. 

As a final example, Figure 10 shows some partial results (10 reconstructed profile curves) from 
an image sequence of a coffee mug. This example demonstrates that our method can handle objects 
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Figure 6: Sample real image sequences used for experiments 
(a) dodecahedron (b) soda can (c) coffee (d) tea. 
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Figure 8: 3D reconstructed points 
(a) dodecahedron, (b) soda can, (c) coffee, and (d) tea. 
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c d 
Figure 9: Profile and epipolar curves 

(a-b) soda can (c-d) coffee. 
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algorithm 


n 




n f >3 


n f >l 


n f > 7 A a 2 x < 0.5 


smoother 


7 


0.0 


.0074 (77%) 


.0046 (45%) 


.0044 (38%) 


smoother 


7 


0.1 


.0114(74%) 


.0054 (41%) 


.0051 (36%) 


batch 


7 


0.0 


.0042 (79%) 


.0036 (56%) 


.0035 (43%) 


batch 


7 


0.1 


.0074 (77%) 


.0054 (53%) 


.0051 (42%) 


batch 


3 


0.0 


.0008 (77%) 






batch 


3 


0.1 


.0159 (75%) 







Table 1 : Root median square error and percentage of edges reconstructed for different algorithms, 
window sizes (n), input image noise cr ; , and criteria for valid estimates (nf-. minimum number of 
frames in fit, o^: covariance in local x estimate). 

These errors are for an ellipse whose major axes are (0.67, 0.4, 0.8) and for a 128 x 120 image. 



with interior holes, since we are not limited to only following the external silhouettes of the objects. 
In future work, we plan to study the events which occur when multiple silhouette curves obscure 
each other in the image sequence (which corresponds to points of bitangency in 3D). 

10 Discussion and Conclusion 

This paper extends previous work on both the reconstruction of smooth surfaces from profiles 
(edge-based multiframe stereo) and on the epipolar analysis on spatiotemporal surfaces. The 
ultimate goal of our work is the construction a complete detailed geometric and topological model 
of a surface from a sequence of views together with an estimate of uncertainty. Towards this end, 
our observations are connected by tracking edges over time as well as linking neighboring edges 
into contours. The information represented at each point includes the position, surface normal, and 
curvatures (currently only in the viewing direction). In addition, error estimates are also computed 
for these quantities. Since the sensed data does not provide a complete picture of the surface, 
e.g., there can be self-occlusion or parts may be missed due to coarse sampling or limitations on 
the camera trajectory, it is necessary to build partial models. In the context of active sensing and 
real-time reactive systems, the reconstruction needs to be incremental as well. 
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Because our equations for the reconstruction algorithm are linear with respect to the measure- 
ments, it is possible to apply statistical linear smoothing techniques, as we have demonstrated. 
This satisfies the requirement for incremental modeling, and provides the error estimates which are 
needed for integration with other sensory data, both visual and tactile. The application of statistical 
methods has the advantage of providing a sound theoretical basis for sensor integration and for the 
reconstruction process in general [Szeliski, 1989; Clark and Yuille, 1990]. 

In future work, we intend to develop a more complete and detailed surface model by combining 
our technique with regularization-based curve and surface models. We also plan to investigate 
the integration of our edge-based multiframe reconstruction technique with other visual and tactile 
techniques for shape recovery. 
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